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We examine the Gross-Pitaevskii (GP) model for Bose-Einstein condensates in parabolic traps 
with attractive interactions. The decay of metastable condensates is investigated by applying the 
instanton formalism to the GP field theory. Employing various dynamical trial states, we derive 
within a coherent state path integral approach a collective coordinate description in terms of the 
condensate radius, in agreement with (and extending) earlier results. We then solve numerically 
for the complete instanton field configuration and compare with the collective coordinate approach. 
Adjusting only the effective mass of the collective coordinate, the two schemes are then in good 
agreement. 

PACS numbers: 



I. INTRODUCTION 

It is well known that a Bose gas with purely attrac- 
tive interactions is unstable and will collapse to a state 
with no thermodynamic limit. However, in the pres- 
ence of a harmonic confining potential, Ruprccht et al. 
jl| have shown that a metastable condensate exists, pro- 
vided the particle number N is sufficiently small. The 
condition for a metastable condensate is N < N c , with 
N c oc £/|a|, where £ is the quantum confinement length 
associated with the trap, and a is the s-wave scattering 
length; a < corresponds to an attractive s-wave pseu- 
dopotential. Such is the situation with 7 Li, for which 
a = — (14.5 ±0.4) a B , where a B is the Bohr radius Q. In- 
deed, experiments by Bradley et al. 0] demonstrated the 
existence of metastable 7 Li condensates of limited par- 
ticle number (N < N c ). Experiment and theory are in 
rough agreement on the value of N c . 

The metastable condensate has a finite lifetime - it 
tunnels quantum mechanically to an unstable higher den- 
sity state, which then collapses. The tunneling process 
was described by Stoof || in terms of a collective coor- 
dinate q(t) which parameterizes a Gaussian condensate 
density, 

The quantum mechanics of this collective coordinate 
arises from the more microscopic Gross-Pitaevskii (GP) 
field theory || , which is itself a simplification inasmuch as 
interatomic potentials are therein replaced by an effective 
s-wave pseudopotential g S(r). To compute the tunneling 
rate, one must examine the Euclidean action 5 E for the 
GP model, 
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FIG. 1. The potential energy V(q) derived from the Gaus- 
sian condensate profile of equation hi for three values of v. 
When v < v c — I6/I5 5//4 ~ 0.542 a metastable condensate 
exists. 

s E /h= J dT J d3r j^^ + ^W- VV> 

+ I r 2 ^V+is(#) 2 }, (2) 

Here and below distance is measured in units of I = 
\fUjrwjj, where m is the boson mass and u is the trap 
frequency, time in units of and the fields ip and ip in 
units of £ -3 / 2 . The coupling constant is then g — ina/i, 
where a is the (negative) s-wave scattering length. Pe- 
riodic boundary conditions are enforced over the imagi- 
nary time domain r £ [— ^ (3, where (3 is the inverse 
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temperature (units of TkS), Writing the complex order 
parameter field ip = yfp exp(i#) in terms of its ampli- 
tude and phase, Stoof formally integrates out the phase 
field 9(r, r), resulting in an effective action S E [p(r, t)] in 
terms of the density alone. Using the trial density from 
equation |], the action functional becomes 



Gaussian Trial State ; i/=0.40 



S E /h = N dr 



-to q 



V(q) 



where m* is an effective mass and 

3 1 3 2 v 1 



2tt q 6 



(3) 



(4) 



and v — N\a\/i = N\g\/An. A large value of TV justifies 
the semiclassical treatment of this problem. The above 
effective potential V(q) reflects the kinetic energy of con- 
finement, the trap potential, and the attractive interpar- 
ticle potential, respectively. Provided v < v c = 16/15 5 / 4 , 
V{q) has a local minimum at q = q , corresponding to 
a metastable condensate. This state is isoenergetic with 
a denser condensate at q = q n < q , and by considering 
motion in the inverted potential —V(q) one constructs 
the usual bounce trajectory q h (r) [Q and from it the 
tunneling rate T oc exp(— AS E /h), where 



AS E = S E [q b (r)] - S E [q a ] 



NhVS^ f % dq \jv(q)-V{q ) 



(5) 



In this paper we shall use the formalism of coherent 
state path integration [|| to treat the problem of the col- 
lapsing metastable condensate. The instanton configura- 
tion for the GP field theory of equation || is described by 
two complex fields ip(r,r) and ty. Restricting our atten- 
tion to trial field configurations, we recover the collective 
coordinate effective action of equation |^. 



II. DYNAMICS AND QUANTUM TUNNELING: 
COLLECTIVE COORDINATE APPROACH 



Our goal is to numerically solve for the instanton of the 
GP action in the zero temperature limit by numerically 
solving for the saddle point of the action-extremizing 
equations. First, though, we concern ourselves with trial 
functions tp(r,T) which yield a one-parameter effective 
action in the form of equation ^|. In the zero tempera- 
ture (/3 — > oo) limit, the initial and final configurations, 
■0(±i/3, r) should represent a metastable condensate. At 
the temporal midpoint r — 0, ^>(0, r) represents the state 
of the field after it emerges from the tunnel barrier, i.e. 
a denser condensate whose real time evolution (accord- 
ing to the nonlinear Schrodinger equation) is unstable 
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FIG. 2. Contour plots of H(p,q) (top) and H(ip,q) (bot- 
tom) for the Gaussian trial states (y = 0.40). Solid lines 
represent highs and dotted lines lows. (In the top graph, the 
low contours increase linearly while the high ones increase 
quadratically. In the bottom graph, the situation is reversed.) 



towards collapse. To interpolate between these states, 
let us make the dynamical Ansatz 



ip(r, t) — A exp(— a r 2 ) 
t/> = A cxp(— ar 2 ) 



(6) 
(7) 



where A, A, a, and a are functions of the time. Number 
conservation relates these quantities, since 



N = d 3 ripip = 



tt 3 / 2 AA 
(a + a) 3 / 2 

Inserting our Ansatz into equation ||, we find 
S E /h = N i dr 



(8) 



2 a + a 

3 1 



a + a 
v 



4 a + a 



(a + a) 3 / 2 



(9) 



plus a term N Jdr d T In A which vanishes owing to peri- 
odicity. 

We emphasize that during the instanton event, a and 
a are not related by complex conjugation. This feature 
of the instanton solution is familiar, for if one treats the 
problem of quantum tunneling in a potential V(x) using 
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FIG. 3. Numerically obtained results for the energy V(q ) 
of the metastable state and action AS E /Nh as a function of 
v. Results are shown for the Gaussian and the Exponential 
trial condensate functions. 



coherent states, the saddle point path is determined by 
the equations 



dz/dr — —dH/dz 
dz/dr = +dH/dz , 



(10) 
(11) 



where z — (x + ip)/\/2, z — (x — ip)/>/2, and H 



V 



V(x). During the instanton event, z ^ z* - indeed 



z = z*(— r) - so the momentum p(r) is imaginary. 
In fact, we can change variables to (q,p), where 



(12) 
(13) 



so that 



S E /h = N Jdr^- 



i{pq-qp) + -{p 2 + q 2 ) 



+ 4« - 



/2vr 



Let us define the Hamiltonian 



3 

H(p, q) = -p 2 + -q 



3- 2| 3 



(14) 



(15) 



The equations of motion are then 



from which we see that ip is real during the instanton 
event, as usual. 

In figure ^| we show contour plots of H(p, q) and 
H(ip 1 q) for v = 0.40. The minimum of H (p, q) occurs at 
(0,(? ) while q 1 is the location of the maximum of V(q). 
The saddle point is traversed during the instanton by p 
becoming imaginary and ip moving around the contour 
line H(ip, q) = V(q ), which encircles the local maximum 
(0,^) of H(ip,q). At t = 0, i/} describes the nucleated 
state; this is of course denser, since q(r = 0) = q n < q . 

Since the dynamics preserves H(ip,q), the instanton 
action is given by 



S E /Nh = 0V(q o ) + -i ipdq 



H(ip,q)=V(q ) 



(18) 



The difference AS E /Nh = S E /Nh - 0V(q o ) is propor- 
tional to the area enclosed by the contour. 

Integrating out the momentum p(r), we obtain the 
action of equations || and [|, with m* = |. Stoof j^] 
states that this is in fact the exact value for Gaussian 
trial states. Within his approximation scheme, he ob- 
tains m* ~ 0.27. 

Similar results are obtained using exponential trial 
functions IBB 



ij)(r,T) = A exp(— ar) 
ip = A exp(— a r) 



from which one derives 

S-pITi — N fdr 1 + -a« + 6(a + Q)" 

E ./ I a + a 2 



(19) 
(20) 



(21) 



With a = | and a = |, one obtains results similar to 
those for the Gaussian model. Integrating out p again 
generates equation || but with m* — 9 and 



V(q)= 1 -\ + 6q 2 -^-± 
8 q z 32 q A 



(22) 



Numerical results (see figure g) show that for any given 
value of v, that the metastable state energy V(q ) is lower 
for the Gaussian model. The computed action AS E /h 
differs little, however (although the maximum value f max 
is model-dependent). 

If q is the condensate radius, the functional form 



- i q = OH/ dp 
3 

- ip= dH/dq , 



(16) V(q) =bq- 2 + cq 2 -dvq- 3 

(17) 

is generic. Setting dV/ dq = gives q = q Q through 
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FIG. 4. Trial state results for q Q (v) (solid) and q n (v) 
(dashed) . 
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chemical potential jj. 

FIG. 5. Scaled particle number versus chemical poten- 
tial for both Gaussian and exponential models. For fi < p c , 
the condensate is unstable; the curves terminate in a dot at 
these points. The numerical results obtained by solving the 
Bogoliubov equations are shown for comparison. Open circles 
denote unstable condensates. 



which is positive for < q < Q = yjbjc. Stable 
condensates have q > q c where f'(q c ) = 0; this gives 
q c = Qj\fh = {fbjhc. Thus v is restricted to lie in the 
range [0, v c ], where 



^ = f{q c ) = T5d{rc 



1/4 



(25) 



Results for the exponential and Gaussian trial conden- 
sates are summarized in the following table: 



Thus, in the vicinity of the critical point q — q c , one has 
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v c 1 288 
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II + 



(29) 



III. BEYOND THE COLLECTIVE COORDINATE 
APPROACH 
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The chemical potential is given by 
1 , -2 7 2 

hence we find 

u = 0: q= yfbfc , fi = 2 Vb~c 
v = v c : q= ^Jbjhc , ^ = M c 



3V5 



Vb~c . 



(26) 

(27) 
(28) 



We now discuss a full solution of the GP equations, 
where all degrees of the condensate are retained. We 
begin with the condensate itself and the real time Bo- 
goliubov equations for excited states. We begin with the 
nonlinear Schrodinger equation for <fi = y/\g\ip, 

id t 4> = - l -v 2 4> + \r 2 ^-(U)^-n4> , (30) 

where \x is the dimcnsionlcss chemical potential in units of 
Tiuj (we now work in the grand canonical ensemble). We 
are interested in stationary states, and we consider only 
real, radially symmetric solutions, for which we write 
4>{r) = R(r)/r. The radial wave equation for R(r) is 



ld 2 R 1 2 

1- (-r 

2dr 2 { 2 



R 3 

M )i?--=0 



(31) 



subject to the boundary conditions R(0) = R(oo) = 0. 
We solved this equation on a one-dimensional grid and 
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FIG. 6. The three lowest values of Q 2 , eigenfrequencies 
of the Bogoliubov equations, as a function of the chemical 
potential. A negative f2 2 indicates an unstable condensate. 



compute the excitation energies; the spectrum was inves- 
tigated in |Q . In figure || we plot the three lowest values 
of SI 2 versus /i. Consistent with previous work we 
find that for /i < 0.4 there is a purely imaginary eigenfrc- 
quency reflecting the instability of the stationary state. 
The inverse of the imaginary part of Q gives the typical 
time for the instability to grow. For larger values of the 
chemical potential (lower values of N) the unstable mode 
crosses the zero mode and the spectrum becomes that of 
a stable condensate, as in the repulsive case. We find the 
stability boundary lies at v = 0.57, in agreement with 
the results of ref . 0] . 



used a relaxation method jfl]]. The size of the system 
was chosen large enough to capture the structure of the 
wave function. We used as many as 1, 000 grid points in 
the range r £ [0, 10]. The results were checked for inde- 
pendence on the system size. This was done for various 
positive values of the chemical potential /i; we obtained 
the ground and several excited states. The number of 
atoms is given by N = 4atv/\g\, with 



/>oo 

v= I drR 2 (r) 
Jo 



(32) 



The dependence of v on the chemical potential is shown 
in figure [3] along with the predictions of the trial conden- 
sate models. Stationary condensates exist for all values 
of fi < 1.5. The stability of these states can only be 
checked by a numerical integration of the time depen- 
dent nonlinear Schrodinger equation (NLSE) 0], or, as 
we did, by computing the spectrum of fluctuations above 
these condensates. 

Once the ground state solution R is found, we solve 
for the spectrum of fluctuations, writing <p = <p + rj, with 
R (r)/r. This leads to the Bogoliubov equations 

L2LHa| 



n u 

-Qv 



hj 2 + 1 -r 2 -^-2\^\ 2 )v-<t>l 2 {r)u CW) 



where 



V(r,t) = u(r) exp(-int) + v* (r) exp(zfi* i) . (34) 



The condensate function </> is taken to be real. Note that 
(it, v) — (4> , — <^> ) is a zero mode. 

Taking advantage of the radial symmetry of </> ), one 
can classify the excited states by angular momentum and 




FIG. 7. Plot of the radial wave function of the bounce for 
fj, = 0.7. The system size was r £ [—12, 12] and r £ [0,3.5]. 
The parameter S (see equation ^) is 0.028. 

We now compute the tunneling rate using the famil- 
iar instanton formalism [Q. In order to compute the in- 
stanton action, we must solve the following equations of 
motion derived from equation ^. Adding a chemical po- 
tential term, we obtain the equations 



a r = -iv 2 0+(ir 2 -^)0 



(35) 
(36) 



subject to periodic boundary conditions on the interval 
t £ [— h(3, were (3 is the inverse temperature (we are 
interested in the /3 — > oo limit). Note that these equa- 
tions yield the imaginary time version of the continuity 
equation, 



with 



id T p + V • j = 



3 = -(W-W) 
2,1 



(37) 



(38) 
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FIG. 8. The metastable state and the nucleated state for 
the bounce of figure |^. The inset shows — ln<^ vs. r 2 ; devia- 
tions from pure Gaussian behavior are apparent. Thin solid 
lines show the corresponding results for periodic solutions. 

Thus, the total particle number N — j<Pr p{r) is con- 
served by the evolution equations, since j vanishes as 
r — * oo. It should be emphasized that equations [55] and 
|36|are not complex conjugate pairs. However, the instan- 
ton solutions of interest all have the symmetry 



4> = <t>*{r,-r) 



where (*) denotes complex conjugation |l4l - Substitut- 
ing this into equation |3^ leads to a nonlocal equation in 
(imaginary) time for (j). Nevertheless, a straightforward 
application of Newton's method may be used to obtain a 
numerical solution || . The equations admit a static solu- 
tion, <fi(r,T) = </> (r), describing the metastable state. In 
the limit (3 — > oo, the nontrivial bounce solution satisfies 

<t>b{ r i T ) — 4>o ( r ) at T = ±2^' wnue describing a denser 
nucleated condensate at r = 0. The tunneling rate, ob- 
tained by summing over all multibounce solutions pi, is 



(40) 



where A 



T = Acxp(-AS E /7i) 
det~ 1 / 2 \S 2 S w ] is the fluctuation determi 



nant prefactor, and AS E is the difference between the 
Euclidean actions of the bounce and static solutions to 
([35|,[56|). Since the energy is preserved by the dynamics, 
we obtain 



A5 E 
Nh 



fd^< 



(41) 



The numerical method used to solve for 0b (r, r) is de- 
scribed in ref. Briefly, we look for real, radially sym- 
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FIG. 9. A comparison of the numerically obtained instan- 
ton action with various trial state predictions of AS E /Nh. 
Filled squares represent numerical data with aperiodic bound- 
ary conditions; open squares denote periodic boundary con- 
ditions. 



metric solutions, <t>b(r) — Rh{r)/r. We discretize space- 
time and write ( |3q , |36| ) as a set of nonlinear equations 
which we solve using the Newton-Powell method Jl7| . 
This process demands an initial guess for the bounce so- 
lution, which we took to be an interpolation between the 



(39) 

static o (r) at 



and a somewhat denser Gaus- 



sian profile at r = 0. It is necessary to choose this guess 
to be asymmetric in imaginary time about r = 0. 



When (3 is finite, periodic solutions will deviate from 
the static 4> at r = ±^/3. We enforced 4>(r 7 —\(3) = 4> (r) 
which results in non-periodicity of the bounce solution. 
We confirm that as we increase (3 that 4>(r, +hj3) ap- 
proaches </>o( r ) while the central structure around t w 
does not change appreciably. It is the central part of the 
bounce which contributes most to (41 ). A measure of the 
aperiodicity is 



5 = 



fdr\R h (r,-±p)-R b (r,+±l3)\ 

JdrR h (r,~\(3) 
o 



(42) 



We increased (3 until 5 < 0.05, at which point we deemed 
our solution a good approximation of the (3 = oo bounce. 
The radial coordinate was chosen to extend from r = to 
a maximum r = L, where L was chosen to accommodate 
the metastable state. We required Rb(L, t) < 0.01 in all 
cases, where typically we took L = 5. 

Another approach is to require periodicity for all /3, 
as is necessary when computing the finite temperature 



G 



instanton solution. In this case, the wavefunctions at r = 
±i/3 agree with each other, but not with the metastable 
solution of the NLSE. Of course agreement is recovered 
in the (3 — * oo limit. In figures ^ and |^ we compare the 
results of these two approaches for (3 — 24. We found the 
aperiodic solution systematically yielded a lower value of 
AS E /Nh. The agreement between the two schemes for 
j3 = 24 is at the 30% level, with better agreement for 
larger values of \x. 

In figure ^ we plotted the radial wave function of 
the bounce, i?b = rtj)^, for = 0.7. We stress that 
0b is interpretable as a condensate wavefunction only 
at t = ±i/3 and at r = 0, where it represents the 
metastable and nucleated condensates, respectively. In 
figure H we plot, also for [i — 0.7, the radial wavefunctions 
i?o(r, — and Rb(r, 0) corresponding to the metastable 
and nucleated states, respectively. Real time evolution 
of the nucleated state using the nonlinear Schrodinger 
equation describes an ever collapsing condensate. This 
is analogous to a tunneling particle rolling down the hill 
after it emerges from the tunnel barrier. In figure [l^ 
we show the radial condensate density, r 2 |0(r)| 2 , devel- 
oping towards collapse as the nucleated state evolves in 
real time according to the NLSE jl5) . 
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FIG. 10. Real time evolution of the nucleated condensate 
(dark solid curve) for t = 0.05, t = 0.10, t = 0.15, t = 0.20, 
t — 0.25, and t = 0.30. N is accurately conserved by the 
integration. The inset shows the collapse of the condensate 
as measured by (r{t)). 

We worked in the range 0.6 < /i < 1.2. We could not 
reach larger values of \i because these demanded a very 
large f3 which increased the calculation time . For the 
smaller values of /i the bounce becomes flatter and loses 
its temporal structure. Eventually AS E /Nh — > as /i 
decreases to its minimum allowed value of ^ min — 0.4. In 



this region, we only obtained convergence to the trivial 
(metastable) solution. 

A comparison of our numerical results with those of the 
various collective coordinate theories is shown in figure |^. 
Apparently Stoof's approximation to* = 0.27 is in much 
better agreement with the data than to* = I derived 
from the Gaussian trial states of equations (y, []]). No 
other parameters have been adjusted. 



IV. CONCLUSIONS 

In this paper, we investigated the decay of metastable 
Bose-Einstein condensates with attractive interactions 
confined in a parabolic trap. The decay rate was cal- 
culated using the instanton formalism, applied to the 
Gross-Pitaevskii action, 5 E [^, The field theory can 
be reduced to an effective quantum mechanics problem 
by focusing on a collective coordinate which describes 
the width of the condensate, as was first done by Stoof 
f|. We went beyond the collective coordinate approach, 
numerically solving for the instanton field configuration 
ip(r, t). This was accomplished by rendering the nonlin- 
ear partial differential equations for tf> and ^ as a large set 
of nonlinear equations on a space-time lattice and solv- 
ing them via the Newton-Powell method. We can thereby 
derive the nucleated state wavefunction, which describes 
the condensate at the instant it emerges from the tun- 
neling barrier. The nucleated state, evolved according to 
the real time nonlinear Schrodinger equation, collapses 
to an arbitrarily dense state for which the attractive GP 
model is no longer a good approximation. While the nu- 
merically obtained metastable and nucleated states show 
deviations from pure Gaussian behavior, a comparison 
of the leading contribution to the semiclassical tunnel- 
ing rate reveals moderate agreement between numerical 
values and those obtained from Gaussian approximations 
to the (dynamical) condensate wavefunction. Adjusting 
only the collective coordinate effective mass to* puts the 
two schemes in quite good agreement. Thus we conclude 
that the collective coordinate description adequately cap- 
tures the essential physics of the tunneling process. 
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